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Abstract 


This  paper  extends  earlier  results  in  the  area  of  variance  reduc¬ 
tion  techniques  applied  to  simulation  on  a  computer.  In  particular,  it 
views  the  antithetic  sampling  technique  as  a  combination  of  rotation 
and  reflection  sampling  on  a  circle.  The  covariance  structures  induced 
by  the  techniques  separately  and  together  are  derived  and  conditions 
under  which  they  are  optimal  sampling  plans  are  described.  Rates  of 
convergence  for  the  variance  of  the  sample  mean  are  given  for  bounded, 
continuous  and  discrete  random  variables  and  for  unbounded  continuous 
random  variables  with  special,  although  commonly  encountered,  structure. 

The  advantage  of  reflection  (basic  antithetic)  sampling  is  greatest 
when  a  certain  symmetry  property  holds.  Rotation-reflection  sampling 
is  superior  to  rotation  sampling  alone  for  continuous  functions.  In  the 
bounded  continuous  case,  convergence  is  faster  with  rotation-reflection 
sampling.  In  the  unbounded  continuous  case,  the  results  show  that 
rotation-reflection  sampling  speeds  convergence  to  the  large  sample 
convergence  rate  achievable  with  rotation  sampling  alone.  For  the  dis¬ 
crete  case,  rotation  sampling  does  as  well  with  regard  to  convergence 
as  rotation-reflection  sampling  does.  However,  analysis  of  the  discrete 
case  shows  that  a  sample  size  n  may  be  considerably  better  than  another 
n'  although  n*  >  n  . 
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1.  Introduction 


Among  statistical  topics  that  arise  in  computer-based  simulation 
experimentation,  variance  reduction  has  long  occupied  a  central  position, 
conceptually  if  not  in  practice.  Variance  reduction  denotes  the  objec¬ 
tive  of  adding  procedures  to  an  experiment  that  allow  one  to  obtain  a 
specified  accuracy  for  less  cost  than  one  can  achieve  in  their  absence. 
Conversely,  for  a  specified  cost,  a  variance  reduction  technique  enables 
one  to  estimate  parameters  more  accurately  than  one  can  without  such  a 
technique.  The  subject  is  not  new,  topical  publications  having  appeared 
over  twenty  years  ago  (Hammersley  and  Morton  (1956),  Hammersley  and 
Mauldon  (1956),  and  Handscomb  (1958)).  Most  textbooks  on  simulation 
acknowledge  the  relevance  of  the  issue  { e.g Emshoff  and  Sisson  (1970), 
Fishman  (1973),  Fishman  (1978),  Gordon  (1969)  and  Naylor  et  al.  (1965)). 

Unfortunately,  attempts  to  implement  this  noble  concept  in  practice 
have  produced  few  documented  cases  of  success.  One  notable  exception  is 
Carter  and  Ignall  (1975).  No  doubt,  a  principal  reason  for  the  paucity 
of  success  arises  from  the  limited  development  of  the  variance  reduction 
technique  that  has  appeared  in  scholarly  journals  beyond  the  original 
conceptualization  of  Hammersley  and  colleagues.  However,  evidence  of 
change  is  in  the  air.  Recently  Lavenberg,  Moeller  and  Sauer  (1979) 
have  attempted  to  broaden  and  deepen  this  development  with  regard  to  the 
variance  reduction  technique  known  as  the  control  variate  method  as  it 
applies  to  discrete  event  simulation.  They  enumerate  the  do's  and 
dont's  of  the  method  with  examples  that  should  prove  helpful  to  potential 
users.  Also,  Schruben  and  Margolin  (1978)  describe  random  number  stream 
manipulation  techniques  designed  to  induce  variance  reduction. 
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The  purpose  of  the  present  paper  is  to  extend  the  development  of 
the  antithetic  variate  method  of  variance  reduction,  a  procedure  first 
described  in  Hammersley  and  Morton  (1956).  Our  results  considerably 
augment  those  of  earlier  work  in  this  area  and  were  motivated  by  obser¬ 
vations  made  in  Fishman  (1979),  which  described  an  application  of  anti¬ 
thetic  variates  to  population  growth  simulations.  Two  examples  illustrate 
conceptually  the  value  of  the  method  of  antithetic  variates.  Firstly, 
consider  the  evaluation  of  the  integral 

1 

4>  =  /  g(x)  dx 
0 


where 


g2(x) 


dx  < 


If  an  analytical  solution  is  unavailable,  one  can  turn  either  to  numerical 
integration  or  to  the  Monte  Carlo  method.  Let  U-j  • . . .  ,U  denote  a 
sequence  of  independent  observations  from  the  uniform  distribution  on 
[0,1]  .  Let  U(0,1)  denote  this  distribution.  Then  an  unbiased  estimator 
of  4>  is 

.  1  n 
$  =  -  Y  g(U . ) 
n  n  L i  j' 

J 


with 


var  $n  =  0(l/n) 

so  that  the  standard  error  of  $n  decreases  as  0(1 /n^ )  . 

This  random  sampling  is  the  most  elementary  application  of  the  Monte 


-r 


rr 
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Carlo  method.  Variance  reduction  techniques  denote  the  use  of  more 
advanced  sampling  designs  intended  to  speed  the  convergence  of  var  $>n  . 
In  particular,  the  method  of  antithetic  variates  aims  at  inducing  a 
joint  distribution  among  for  which 

var  =  o(l/n) 

while  preserving  the  marginal  distributions  as  U(0,1)  ,  which  guaran¬ 
tees  the  unbiasedness  of  . 

n 

As  a  second  example,  one  may  wish  to  apply  variance  reduction 
techniques  to  a  discrete  event  simulation.  Consider  the  single  server 
queue  with  i.i.d.  interarrival  times  A^.A^,...  ,  i.i.d.  service  times 
SpS^,...  ,  {A^ }  and  {S^  independent,  and  mean  waiting  time  y  . 

From  Findley's  equation  one  has  for  the  waiting  time  of  the  ith 
completion 


Wi  =  max(0,  W._1  +  S.  -  A.. ) 


i  =  1 ,. . .  ,m 


on  a  run  terminated  after  m  completions.  Also,  let 


W 


(J  =max(0,  ♦  S(J  -  Ay) 


j  =  1 .  ,n 


denote  the  waiting  time  of  completion  i  on  the  jth  of  n  replica¬ 
tions.  As  an  estimator  of  y  one  has 

n 


Vn  "  n  W*j 


W  .  _  i  l  w..  . 

•j  =  rn  z'i_1 


where 
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If  replications  are  independent,  then  var  n  =  0(l/mn)  .  One 
possible  application  of  the  method  of  antithetic  variates  might  be  to 
induce  joint  distributions  among  A.  , ,  ...»  A.  and  among  S.  ,,  .... 

l  5  i  i  »ii  i  >  i 

S.  so  that 
l  ,n 

var  um  n  -  0(l/m)  o(l/n)  . 

Here  entries  within  a  column  of  the  array 


1,1 

A2,l 

A  i 
m,l 

1,2 

A2,2 

Am,2 

• 

• 

• 

1  ,n 

A«  n 
2,n 

. . . 

Am,n 

are  correlated,  but  entries  within  a  row  are  independent.  A  similar 
characterization  applies  to  service  times. 

At  this  point,  it  is  important  to  recognize  that  the  direct  appli¬ 
cation  of  the  antithetic  method  to  be  described  here  does  not  necessarily 

achieve  o(l/n)  in  var  y  ,  a  well-known  fact  in  multivariate  Monte 

m,n 

Carlo  sampling.  For  example,  see  Hammers  ley  and  Handscomb  (1964). 
However,  a  comprehensive  understanding  of  how  the  technique  works  in 
univariate  problems  is  a  prerequisite  to  devising  methods  that  will 
achieve  the  desired  effect  in  multivariate  problems  such  as  the  aforemen¬ 
tioned  queueing  simulation. 

The  formal  concept  of  the  antithetic  variate  method  first  appeared 
in  Hammers ley  and  Morton  (1956).  Two  subsequent  papers.  Hammers ley  and 
Mauldon  (1956)  and  Handscomb  (1958),  demonstrated  a  certain  optimal 


/  ' 
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property  of  the  method.  Andreasson  (1972)  and  Andreasson  and  Dahlquist 
(1972)  introduced  the  formalisms  of  group  representation  as  a  way  of 
analyzing  potential  antithetic  sampling  designs,  and  Roach  (1973) 
attempted  to  formalize  the  topic  as  a  transportation  assignment 
problem.  The  present  paper  examines  and  extends  the  formulations  in 
these  early  papers  into  an  account  that  sheds  considerable  new  light  on 
the  antithetic  method  and  how  it  works.  Section  2  reviews  the  formalisms 
of  the  antithetic  variate  method.  Section  3  describes  a  procedures, 
based  on  rotation ,  for  collecting  n  >  2  antithetic  replications  that 
lead  to  considerably  greater  accuracy  per  unit  cost  than  the  tradition¬ 
ally  recommended  antithetic  variate  method  for  n  =  2  allows.  It  also 
describes  several  examples  that  reveal  how  this  rotation  sampling  performs 
in  selected  situations.  Section  4  describes  a  procedure,  based  on 
rotation  and  reflection ,  that  in  certain  cases  improves  on  the  method 
of  Section  3  ,  and  illustrates  its  app1!' cation  to  some  of  the  examples 
in  Section  3  .  Section  5  describes  one  circumstance  in  which  the  results 
derived  here  apply  to  a  single  server  queueing  simulation.  Both  rotation 
and  reflection  sampling  designs  make  clear  the  value  of  continued  study 
of  these  procedures. 

2.  Basic  Antithetic  Sampling 

Consider  the  random  variables  n-|..--»nn  and  suppose  one  forms  the 
quantity  h(nj»...,nn)  and  uses  it  as  an  estimate  of  an  unknown  quantity 
4>  .  If 

n_1  E  h(n1 ,. . .  ,nn)  =  <j> 

the  estimator  is  unbiased.  Moreover,  a  low  value  for  var  h(n-j ,. • . ,nn) 
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o 

relative  to  <p  indicates  high  reliability  for  h(n-) ,  •  •  •  ,nn)  as  an 
estimator  of  4>  .  An  important  subclass  of  interest  is  the  separable 
function 


h  (n-j ,  -  -  -  ,nn )  =  h-j  (n-j )  +  h2(n2)  +  •••  +  hn(nn)  • 


Gi  ven 


h,,...,h  ,  the  marginal  distributions  of  and  the 

In  in 


condition  E  /  h(n-)  =  n<<>  ,  one  can  concentrate  on  choosing  a  joint 

j  =  l  J 

distribution  for  n-|»...»Tin  to  promote  reliability  without  concern  for 
bias.  Working  in  a  different,  but  related,  function  space  facilitates 
this  choice. 

Let  ni  have  the  cumulative  distribution  function  (c.d.f.)  F. 
with  inverse  distribution  function 


G j ( x )  =  inf[y:  F^(y)  >  x,  0  <  x  s  1]  . 


Let  U,  U-j ,  ...,  Un  denote  uniform  deviates  and  defi 


ne 


*  yvuj)]  ■  w 


Then  the  estimator  of  $  of  interest  is 


Tn=l  L,,  9J(Uj>  '  (1) 

One  can  now  restate  the  variance  reduction  problem:  Given 

with  E  g.(U.)  =  0  ,  choose  the  joint  distribution  of  U-, , . . .  ,U  to 
J  J  in 

minimize  var  T  . 

n 

At  this  point  the  Antithetic  Variate  Theorem  becomes  salient. 


Theorem  1  .  Define  Q  as  the  set  of  all  functions  for  which 


i.  oj( z )  is  a  1-1  mapping  of  (0,1)  onto  itself, 
ii.  Except  at  a  finite  number  of  points  z  ,  dw/dz  =  1  . 

Also,  I  =  infimum  var  T  over  all  possible  stochastic  and  functional 
n  n 

dependences  among  U-|,...,Un  .  Then 

inf  var[  j-  l  gJwJU))]  =  I  .  (2) 

UjCil  n  j=l  J  J  n 

j  =  1 , . . .  ,n 

For  bounded  g.|,...,g  Hammersley  and  Mauldon  (1956)  give  the  proof 
for  n  =  2  and  Handscomb  (1956)  gives  the  proof  for  n  >  2  .  Recently 
Wilson  (1979)  has  extended  the  theorem  to  unbounded  g^,...,g  . 

Theorem  1  has  profound  implications.  It  says  that  one  can  achieve 
the  infimum  I  by  generating  a  uniform  deviate  l)  and  applying  measure 

preserving  transformations  on  (0,1)  .  As  an  example,  consider  the  case 
of  n  =  2  ,  h-j(x)  =  h2(x)  and  monotone,  g^(x)  =  g2(x)  ,  and  G2(y)  = 
G-|(l  -  y)  •  Then  the  sampling  design  uj^(U)  =  u>2(U)  =  U  gives 


T2  =  \  C^i  (G] (U))  +  h2(G2(U))] 


(3) 


=  \  Cg-j (u )  +  g-j  ( l  -  u)] 

for  which  var  T2  =  I2  .  In  the  case  h^(x)  =  x  ,  g-|  =  G-j  the  minimal 
variance  implies  that  no  other  method  of  generating  and  n2  produces 
a  more  negative  correlation,  a  result  well  kno^’n  in  probability  theory. 
See  Hoeffding  (1940),  Frechet  (1951),  Mardia  (1976)  and  Whitt  (1976). 

It  is  this  form  of  basic  antithetic  sampling  (u)-j(U)  =  to2  (U))  that 


i 


textbooks  on  simulation  usually  describe. 

The  problem  that  now  arises  is  to  choose  {w-(U);  j  =  l,...,n}  that 

J 

achieves  the  infimum  of  var  Tn  for  n  >  2  .  This  is  not  a  simple  problem 
nor  do  we  claim  to  have  solved  it  entirely.  However,  our  results  are 
encouraging.  Section  3  describes  the  concept  of  rotation  sampling  and 
shows  its  optimality  under  specified  conditions.  Section  4  then  combines 
basic  antithetic  sampling  with  rotation  sampling  into  a  rotation  and 
reflection  sampling  scheme  that  considerably  accelerates  the  convergence 
of  var  Tn  . 

3.  Rotation  Sampling 

The  task  of  selecting  among  alternative  measure  preserving  trans- 
mations  on  [0,1)  can  be  simplified  at  the  outset  by  considering  two  par¬ 
ticular  sets.  Firstly,  consider  transformations  of  the  form 

<4U)  =  U  0  s  U  <  $ 

=  1  +  6  -  U  0  <  U  <  1  0  <  0  <  1 


Figure  la  shows  an  example.  Since  these  fail  to  satisfy  point  ii  of 
Theorem  1  ,  we  omit  them  from  further  consideration.  As  a  second  alter¬ 
native,  consider  the  transformations 


u>(U)  =  u 

=  U  +  63  -  62 
■  u  +  e]  -  e2 
=  u 


0  s  U  <  31 
61  ^  U  <  02 
<63 

63  s  U  <1 
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for  fixed  0  <  0^  s  ^  6  S3  <  1  •  Note  that  <*>  covers  the  unit  interval 
in  non-overlapping  segments.  Figure  lb  shows  an  example.  Observe  that 
mappings  of  this  form  have  several  constants  to  be  evaluated,  thereby 
adding  to  the  selection  problem. 


Figure  la  Figure  lb 


Here  we  confine  our  attention  to  the  set  of  one-parameter  transformations 
«j(U)  >  U  •  ■  U  +  0  <  U  <1  -  0j 

■  u  +  0,  -  1  l  -  ei  <  u  <  1  (4) 

J  w 

0  S  6j 

for  j  ■  l,...,n  .  Since  these  transformations  constitute  rotations  on  the 
unit  circle,  we  refer  to  (4)  as  rotation  sampling.  For  convenience  of 
exposition,  assume  that  a)  have  cornnon  c.d.f.  F  with 

corresponding  inverse  distribution  function  G  and  b)  h^(x)  =  ...  =  hn ( x ) 

=  h(x)  .  Since  gj  *  h^G^  ,  it  follows  that  g^  =  ...  =  gn  =  g  . 
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1  2 

Lastly,  assume  that  c)  /  g  (u)  du  <  00  . 

0 

One  can  now  write  (1)  as 

T„  -  I  9(U  «  Oj)  (5) 

and  let 

P(6)  =  E  g(U)  g(U  «  0)  -  4>2  .  (6) 

Among  the  properties  that  follow  from  (5)  and  (6)  are 
Property  1  .  (unbiasedness)  E  g(U  ®  6.)  =  E  g(U)  =  0  . 

Property  2  .  (continuity)  The  function  P  is  continuous  on  (-=°,oc)  . 

Property  3  .  (differentiability).  If  g  is  continuous  on  [0,1]  , 
then  P  is  differentiable  on  [0,1]  .  If  g  is  continuous  but 
unbounded  at  u  =  1  ,  then  P  is  differentiable  on  (0,1)  .  If 
n-|,...,nn  have  a  discrete  or  mixed  marginal  distribution,  then  P  has 
nondifferentiable  points  in  (0,1)  . 

Property  4  .  (periodicity)  P(9)  =  P(9  mod  1)  0  £  (-“',«)  . 

Property  5  .  (symmetry  about  0  =  1/2)  P(0)  =  P(1  -  6)  . 

Property  6  .  (symmetry  about  0=0)  P(-e)=P(0). 

1 

Property  7  .  (exhaustiveness  )  /  P(6)  de  =  0  . 

0 

Property  8  .  (stationarity) 

cov[g(U  ®  Q-),  g(U  ®  0k)3  =  E  g(U  9  0j )g(U  ®  6k)  -  <t>2  =  P(0j  -  ek) 

0  s  ej5ek  . 

Property  9  .  (upper  bound)  P (0)  =  P(l)  -  P ( e )  9  £  (0,1)  . 

Property  10  ■  (lower  bound)  If  P  is  convex  on  [0,1]  ,  P(l/2)  <  P(e)  . 
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These  properties,  especially  those  relating  to  stationarity  and  symmetry, 
prove  useful  when  selecting  0^,...,0  to  minimize  var  Tn  . 

This  problem  can  now  be  formulated  as 

n-1  n 


mi 

0 


inimize  V(0,,...,0  )  =  l  l  P(9  -  e  ) 

1 » • . • >6  k=l  j=k+l  J  k 


1  ’  ’  "  ’"n 
subject  to 


o  <  e1 
ej s  Vi 

8n  s  1  • 


j  =  1 


(7) 


.n-1 


Expression  (7)  is  equivalent  to  minimization  of  the  average  correlation 
coefficient  of  h(n-| ) ,. . .  ,h(nn)  .  This  formulation  leads  to  Theorem  2. 

Theorem  2  .  If  C  is  a  convex  function  on  [0,1]  and  symmetric  about 
z  =  j  ,  then  for  given  n  >  2  z*  =  (~ ,  . . . ,  ^-)  is  an  optimal  solution 
of  the  optimization  problem: 


n-1 

n 

min 

w(z)  -  l 

l  C(z.  + 

1 . zn-l 

)  i=l 

j=i+l  1 

n-1 

subject  to 

<  1 

+ 


) 


(8) 


i  =  1 .... .n-1  . 


See  che  Appendix  for  the  proof. 

Letting  C  =  P  tells  us  that  the  assignment 
for  j  =  l,...,n  gives 


*  J i  “  1  * 
6j  =  E  zk 

J  k=l  K 


(j-D/n 


v(er...,en)  *  v(e] .... .en)  . 

★  ★  * 

The  assignment  0n  =  (0j,...,0  )  leads  to  considerable  convenience. 


i 


mmz'-r 


fr 
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★  ★ 

In  particular,  {oj- (U )  =  U  ffl  6,;  j  =  form  a  finite  cyclic 

J  J 

abelian  group.  Define  P|-j_j|  =  p(0-j  _  ej)  •  Then  { 9 (U  ®  9 j ) }  has 
the  covariance  matrix 


(9) 


Here  row  j+1  is  row  j  with  elements  shifted  one  position  to  the  right 
and  the  right-most  entry  in  row  j  assigned  to  the  left-most  position  in 
row  j+1  .  A  matrix  with  this  property  is  called  a  airaulant.  Its  kth 

eigenvector  is  |e''^7r^^n  ;  j  =  0,...,n-l}  where  i  =  /T  ,  which  gives 
the  eigenvalues 


p  2ri(k-l)j/n 
j 


k  =  1 . n  ,  (10) 


In  particular,  note  that  the  unitary  matrix  Vn  =  ||  e^^  ^J//n|| 

orthogonal izes  £n  and  that  =  Tp_jc+2  n  f°r  ^  =  2,...,n  . 

★  , 

Regardless  of  whether  or  not  e  is  optimal,  the  resulting  symmetry 

★ 

in  Pp  affords  an  understanding  of  the  rate  of  convergence  of  var  Tn 
with  n  where 

T*  =  l  l  G(U  6  (j-l)/n)  . 
n  n  j=l 


*T r 
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Theorem  3  . 


If  one  uses  the  transformations  {w.;  j 

J 


n-1 


var 


Tn  '  S'  F.0  P(n>  ■  Tl,n/n  • 


1 , . . .  ,n }  ,  then 


01) 


Proof  .  Observe  that 


var 


1 

n 


n  n 

I  l 

k=l  j=l 


I  k- j 


where  the  summation  is  over  all  elements  in  P  .  Since  each  row  of 

~n 

£n  contains  the  same  elements,  it  follows  that  the  summation  over  all 

elements  in  P  is  equivalent  to 
~n  n 

*  i  n-1  .  n-1  . 

var  T  =  4  (n  I  P.)  =  -  T  P(J-)  .  (12) 

"  „2  j=0  J  "  j=0  " 

From  (10)  with  k  =  1  ,  it  is  clear  that 

varTn'Tl,n/n  ' 


The  convergence  problem  now  becomes  one  of  showing  how  x, 

I  »n 

behaves  as  n  -*■  »  under  alternative  restrictions  on  g  .  To  put  this 
problem  in  perspective,  observe  that 

van  T*  =  J  ll  1  P(^-)  +  X  P(0)  +  X  P(l)]  -  f]  P(0)  d0 
n  n  j=l  n  L  c  0 

★ 

so  that  one  can  interpret  var  Tn  as  the  error  incurred  in  using  the 
trapezoidal  rule  to  approximate  the  integration  of  P  over  [0,1]  . 

Theorem  4  .  If  one  uses  the  transformation  {w.;  j  =  l,...,n}  then 

J 

★ 

Tn  is  the  minimum  variance  unbiased  estimator  of  <j>  for  fixed  n  . 


Proof.  Consider  the  estimator 


Tn  =  l  c.  g(U  e  e.) 
j=l  J  J 


n 

o  .  =  1  . 

j=l  J 

Let  c„  =  (c, , . . . ,c  )  .  In  order  for  T  to  be  the  minimum  variance 
~n  i  n  n 

unbiased  estimator  of  <p  ,  one  needs 


~n 


Q 


1 


where  1^  is  an  nxl  vector  of  ones.  Since  Vp  is  a  unitary  matrix, 

we  have  V”1  =  vj  and  t  =  VT  V„  where  r  is  a  diagonal  matrix 

~n  ~n  ~n  ~n  ~n  ~n  ~n  3 

with  xk  n  in  row  k  and  column  k  .  Then 


C„  *  til  V„  vj  Xh)'1  iT  V0  t;1  vj  =  1  (1 . 1) 


/SJ  ★ 

so  that  Tn  =  ,  which  proves  the  theorem.  This  result  is  also  noted 

in  Andreasson  (1972). 

★ 

An  equivalent  representation  for  Tn  proves  useful. 

★ 

Lemma  5.1  .  For  the  rotation  scheme  0 
- 


★ 


1 

n 


(13) 


where  £  is  from  U(0,1)  . 


5- 


Proof  .  One  has 


n-1 


t  =  -  y  gfu  e 

n  n  L  n-> 


n-1 


[  J=0  n=m+l  j 


(14) 


where  m  =  [ml  -  U)J  .  Let  £  =  1  -  n(l  -  U)  +  m  =  nU  mod  1  so  that 


n-1 


\--Vj  + 1  '  g  {^jr1) 

n  n  [  j=0  n  j=m+l 


n-1  _,i  n-m-2  r+t,. 

I  st^}  +  l  g(^ 

k=n-m-l  n  k=0  n 


(15) 


J  f 1  ■ 

n  k=0  n 


Clearly  £  is  from  U(0,1)  .  Expression  (13)  is  identical  with  the 
Hammersley  and  Morton  (1956)  formulation  for  n  >  2  .  Their  convergence 
results  make  use  of  the  Euler  summation  formula  (see  Fort  1948,  p.  53) 


i  n_1  y  +  1  m 

l  l  (  ■  /  9(t>  «  ♦  l 

n  j=0  n  0  k=l 


Bk(x)[q(k~1)(l)  -  q(k~1)(0)] 
k!  nk 


+  o(l/nm)  (16) 

where  0  <  x  <  1  and  (x )  denotes  the  kth  Bernoulli  polynomial 
for  an  arbitrary  function  q  whose  first  m  derivatives  exist.  Note 
that  (12)  and  (13)  both  are  amenable  to  this  representation  subject 
to  the  existence  of  the  appropriate  derivatives. 

★ 

Bounded  g  .  We  now  explore  the  convergence  of  var  Tn  under  alterna¬ 
tive  restrictions  on  g  .  Theorem  5  relates  to  bounded  continuous  g 
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with  finite  first  derivative  and  Theorem  6  ,  to  piecewise  linear  g  with 
finite  discontinuities. 

Theorem  5  .  (Hammersley  and  Morton  1956).  If  g  €  C^[0,1]  ,  then 

a.  T; .  ♦  -  IZillam aMl  a  0(1/n) 

b.  var  T  =  — U-  [g(l)  -  g(0)]2  +  o(l/n2) 

n  12n^ 

where  V  =  nU  mod  1  is  from  U{0,1)  . 

★ 

Proof  .  Result  a  follows  from  substitution  into  (16).  Since  T^ 
is  unbiased,  result  b  follows  directly. 

To  appreciate  the  significance  of  this  result,  one  needs  a  measure 
of  variance  reduction.  One  suggested  measure  is 

VR,  *,  _  variance  without  variance  reduction  technique 
v  ^~n  variance  with  variance  reduction  technique 

Then  for  rotation  sampling  with  bounded  g 

lim  n_1VR(6*)  =  0(1 ) 

n-K» 

so  that  variance  reduction  is  0(n)  . 

Example  5.1  .  Consider  a  Beta  random  variable  with  c.d.f.  F(x)  =  xa 
0  <  x  <  1  and  0  <  a  s  1  so  that  G(U)  =  U^a  .  Let  g  =  G  .  Then 
T*  =  a  +  1  '  (v  "  T/2)/11  +  o(l/n)  and  var  Tp  =  l/12n2  +  o(l/n2)  . 

Observe  that  for  0  <  a  <  1  the  corresponding  p.d.f.  is  unbounded  at 
x  =  0  .  For  the  bounded  case  (a  >  1 )  ,  Theorem  5  does  not  apply, 
and  one  needs  an  additional  result. 

A 

Corollary  5.1  .  If  g  is  continuous  on  [0,1]  ,  then  var  Tn  =  o(l/n)  . 

Proof  .  If  g  is  continuous,  P  CC^O.l]  .  Using  (16)  with  P 

and  x  =  0  ,  one  has 

*  1  B, (0)[P(1 )  -  P(0)] 

var  T  «  J  P(6)  de  +  — - - - +  o(l/n)  . 

n  0 


-17- 


Since  P(0)  =  P(l)  ,  var  Tn  =  o(l/n)  . 

A  somewhat  stronger  convergence  result  than  Corollary  5.1  is 
also  possible. 

1 

Corollary  5.2  .  If  g  is  continuous  on  [0,1]  and  f  g'(x)  dx  <  ®  , 
*  o  0 

then  var  Tn  =  0(l/n  )  . 

Hal  ton  and  Handscomb  (1957)  make  this  assertion.  Huang  (1980) 
gives  a  proof. 

Example  5.2  .  For  the  Beta  case  with  c.d.f.  F(x)  =  xa  0  <  x  <  1  and 

a  >  1  ,  one  has  g;(u)  =  -  u^a~^  which  is  unbounded  at  u  =  0  .  How- 

,  J  a 

★  O 

ever,  /  g'(u)  du  =  1  <  00  so  that  var  T  =  0(1  /n^ )  . 

0  n 

Theorem  6  .  If  g  is  piecewise  linear  with  finite  discontinuities,  then 

var  Tn  =  0(l/n  )  .  See  the  Appendix  for  the  proof. 

Example  6.1  .  Consider  a  Bernoulli  random  variable  with  inverse  distri¬ 
bution  function 


G(U)  =0  0  <  U  <  1-p 

=  1  1-p  <  U  <  1 

and  let  g  =  G  .  Then  P(e)  =  (p  -  e)+  +  (p  +  e  -  1)+  -  p2  0  <  0  <  1  , 

where  x+  =  max(0,x)  .  Figure  2  shows  P ( 9 )  .  Note  that  P  is  convex 

but  not  differentiable  at  0  »  p,  1-p  .  Also 

var  T*  •  lmL2°<Lm  -  .("P.Jfi.1  ’>)  .  0  np  •  integer 

n 

s  — U-  always, 

4n^ 


so  that  variance  reduction  is  infinite  when  np  =  integer  and  otherwise 


I 


-18- 


0  p  1  -  p  19 


Figure  2  P(e)  for  a  Bernoulli  Random  Variable 

is  0(n)  .  Moreover,  perusal  of  Table  1  in  Fishman  (1979)  leads  to 

the  conjecture:  Given  n  as  the  maximal  permissible  sample  size,  then 

using  only  n*  =  min{j:  jp  mod  1  is  a  minimum,  j  =  n}  leads  to 

var  T**  <  var  T*  for  j  *  1 . n  .  For  the  more  general  discrete  case, 

n  j 

Huang  (1980)  shows  that  if  F  assumes  only  rational  values,  there  exist 

★ 

n  's  for  which  var  =  0  . 

Unbounded  g  .  Here  we  use  results  from  generalized  function  theory. 
Consider  the  function 

q(u)  =  ua(l  -  u)br(u)  a,b  $  0  0  <  u  <  1 

where  q  is  integrable  and  the  first  m  derivatives  of  r  exist.  Then 
Lyness  and  Ninham  (1967)  give  the  extended  Euler-Maclaurin  summation 
formula 


-19- 


1 

n 


1 

/  q(u)  du 

0 


m-1 

+  1 

j=0 


j! 


g(-a-j,x)  , 
na+^+1 


(-1) 


J  *1 


(J) 


(1) 


j! 


nb+j+l 


+  0(l/nm)  0  <  x  <  1  (17) 

where  ijjQ(u)  =  (1  -  u)br(u),  (u)  =  uar(u)  and  e(*»*)  denotes  the 
generalized  Riemann  zeta  function.  We  then  have: 

Theorem  7  .  If  g  =  q  and  r  £C^[0,1]  ,  then  var  =  0(l/n2<^  a>) 
if  a  <  b  and  var  Tn  =  0(l/n2^1+b))  if  a  >  b  . 

Proof  .  The  result  follows  directly  from  (17)  with  m  =  1  ,  x  =  nL)  mod  1 

★ 

and  the  fact  that  Tn  is  unbiased. 

Observe  that  variance  reduction  if  0(n1+2a)  for  a  <  b  and 
and  0(n1+2b)  for  a  >  b  ,  so  that  the  efficiency  of  rotation  sampling 
increases  only  if  a  >  -1/2  and  b  >  -1/2  .  But  this  is  precisely  the 
condition  that  assures  a  finite  variance  for  g  . 

Example  7.1  .  Consider  the  Pareto  distribution  with  c.d.f  F(x)  =  1  -  xc 

1  /c 

c  <  0  and  inverse  distribution  function  G(U)  =  (1  -  U)  '  .  Also  let 

g  =  G  .  This  representation  corresponds  to  a  =  0  and  b  =  1/c  in  (17) 
so  that  var  Tn  =  0(l/n2(^  +  1/,c))  which  requires  c  >  -2  (a  >  -1/2) 
to  achieve  a  variance  reduction. 

Other  types  of  unbounded  variation  are  also  possible.  Consider  the 
representation 


q(u)  =  ua(1  -  u)br(u)  In  u 


0  <  u  <  1  . 
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If  q  is  integrable  and  the  first  m  derivatives  of  r  exist,  Lyness 
and  Ninham  (1967)  give  the  extended  Euler-Maclaurin  summation  formula 


where  a,b  <  0  and  the  coefficients  e.(x)  are  independent  of  n  . 

J 

This  gives  rise  to  Theorem  8  . 

Jheorem_8_.  If  q  =  g  and  r  £  C^OJ]  ,  then  var  T*  = 

0((ln  n/n1+a)2)  if  b  >  a  and  var  T*  =  0(l/n2(1+b))  if  b  <  a  , 

Proof  .  The  result  follows  from  substitution  into  (18)  with  m  =  1 

and  x  =  nil  mod  1  and  from  the  unbiasedness  of  T  . 

n 

1+a  2 

The  term  0((ln  n/n  )  )  calls  for  additional  study.  Observe 
that  lim  .  This  implies  that  for  suffi- 

n-[\(kn)UaiVln  "/I  k2  +a 

ciently  large  n  var  T^  =  0(l/n  )  for  a  <  b  .  Here  the  logarithmic 

singularity  slows  the  convergence  rate  for  moderate  n  but  ultimately 
has  no  limiting  effect. 

Example  8.1  .  let  g (U )  =  G(1  -  U)  =  -In  U  so  that  g(U)  is  an  exponen¬ 
tial  random  variable  with  unit  mean.  Here  a  =  b  =  0  ,  r(u)  =  1 

and  var  Tn  =  0((ln  n/n)  )  .  Again  for  sufficiently  large  n  var  Tn  = 

2 

0(l/n  )  ,  the  rate  achievable  for  bounded  functions. 
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4.  Rotation  and  Reflection  Sampling 


As  mentioned  in  Section  2  ,  the  use  of  U  and  1-U  leads  to 
the  most  negative  correlation  between  two  random  variables  when  h  is 
monotone.  However,  the  exclusive  use  of  rotation  overlooks  the  benefit 


of  this  transformation.  To  investigate  this  alternative  more  thoroughly. 


we  consider  n 

=  2m  replications,  where 

(4)  defines 

u>i(U),... 

,um(U) 

and 

Gj+m(x)  =  Gj(l-x)  =  G(x) 

j  =  1,-.. 

,m 

(19) 

Vm(U)  =  co.(U) 

j  =  1,... 

,m  . 

(20) 

Although  (4) 

and  (20)  induce  identical 

c.d.f . 1 s  on 

n1 ,. . . ,nn 

*  "j 

and  nk+m  have  different  joint  distributions  than  the  corresponding  ones 
for  and  for  j,k  =  l,...,m  .  Also,  (20)  meets  the  require¬ 
ments  of  point  i  of  Theorem  1  . 

Note  that  9j(x)  =  9j+mO-x)  =  9(x)  for  j  =  l,...,m  .  Also  note 

that 


1  -  oj.(U)  =  1  -  (U  ®  0.)  =  (1  -  U)  0  (1  -  0-) 
J  J 


=  1  -  0.  -  U  0  £  U  <  1-0. 

J  J 

=  2  -  0j  -  U  l-0j  <  U  <  1  . 


One  can  now  write  Tn  as 


t  =  —  y  f  (u 

n  m  L  - 


j  =  l 


(21) 


where 


f(U  #  Oj)  =  -j[g(U  0  6j)  +  g(l  -  (U  ®  Q.))]  . 


(22) 
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The  first  term  in  the  summand  of  (22)  provides  for  rotation  on  the 
unit  circle;  the  second  provides  for  reflection  on  the  unit  circle. 
Properties  of  interest  include: 

Property  11  .  (symmetry)  f(x)  =  f(l  -  x)  . 

Property  12  .  (unbiasedness)  E  g(l-e..  ®  -U)  =  E  g(U)  =  q  . 

Property  13  .  (stationarity )  cov[g(l-e„.  ®  -U)g(l-6,.  6  -U)]  =  Pd",  -  ■  ) 

J  K  ^  j 

Property  14  .  (stationarity) 

cov[g(u  ®  0j)gn'0k  ®  -U)]  =  Q(ek  -  6j )  6k  >  e. 

=  Q(1  -  0k  +  3 j )  6k  -  6j  ’ 

where 


=  /  g(u  ®  9 ) g ( 1  -  u)  du  -  e  =  6.  -  6.  >  0  . 

0  J 
1  1 

Property  15  .  (exhaustiveness)  [  Q(e)  de  =  f  Q(1  -  9)  de  =  0  . 

0  0 

Property  16  .  (symmetry)  If  g(u)  +  g(l  -  u)  =  2g(l/2)  for  0  <  u  <  1  , 
Q(e)  =  q(l  -  e)  for  e  e  [o,i]  . 

Property  17  .  (lower  bound)  Q(0)  =  Q(l)  >  -P(l)  . 

Property  18  .  (continuity)  Q  is  continuous  on  [0,1]  . 

To  investigate  the  simultaneous  benefit  of  rotation  and  reflection, 
it  is  convenient  to  study 

1  1  2 
R(e)  =  j  E2P(e)  +  Q(e)  +  Q(i  -  0)]  =  /  f(u)f(u  «  0)  du  -  f 

H  0 


- 


0  e  [0,1]  . 
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Note  that 

Property  19  .  (symmetry)  R ( 9 )  =  R(1  -  0)  . 

Property  20  .  (upper  bound)  R(0)  =  R(l)  -  R(a)  tor  0  £[0,1]  . 

Property  21  .  (differentiability)  If  g  £  C^OJ]  ,  R  (0)  = 

r'O)  =  o  . 

Proof: 

1-0  1  o 

R  (e)  =  /  f(u)  f (u+0)du  +  /  f(u)f(u+e-l)du  -  $ 

0  1-0 

1-9  1 

R'(0)  =  f(l-0)[f(O)  -  f(l)]  +  /  f(u)f ' (u+0)du  -  /  f(u)f'(u+e-l)du 

0  1-0 

1 

R'(0)  =  f(l)[f(0)  -  f (1 )]  +  /  f(u)f'(u)du 

0 

By  property  11,  f(0)  =  f(l)  and 

4  1 

/  f(u)f'(u)du  =  -  /  f(u)f'(u)du 
0  4 

so  that  R'(0)  =  0  .  By  property  19  ,  R'(0)  =  -R'(l)  .  Note  that  this 

result  does  not  necessarily  apply  if  R'(0)  does  not  exist  everywhere 
on  [0,1]  . 

Property  21  prevents  us  from  deriving  a  result  for  R  comparable 

to  Theorem  1  .  Since  R'(0)  =  R'(l/2)  =  r'(1)  =  0  ,  R  is  not  convex. 

★ 

Nevertheless,  the  choice  of  0n  has  highly  beneficial  properties  which 


the  next  several  theorems  describe. 
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let 


1 

m 


m-1 

l  f(U  ©  j/m) 
j=0 


n  =  2m 


so  that 


m-1 


**  1  1 

,dr  Tn  =  m  ^  R(j/m)  ’ 


j=0 


since  it  is  easily  seen  that  |f(U  ©  j/m);  j  =  0,...,m-l}  is  a  cyclic 
stochastic  process. 


Theorem  9 


If  property  16  holds,  then  var  Tn  =  0  . 


Proof  .  Since 

g(U  ©  j/n)  +  g(l  -  (U  ©  j/n))  =  2g(l/2)  j  =  0,...,m-l  . 

T**  =  9(1/2)  and  var  T**  =  0  . 

Bounded  g  . 

Theorem  10  .  If  g€C^[0,l]  then 

a.  Tn  =  <p  +  o(l/n)  . 

b.  var  Tn*  =  o(l/n2)  . 

Proof  .  Let  f(u)  =  q(u)  as  in  f!6)  and  take  x  =  V  =  nU  mod  1  . 

By  property  11  f (0 )  =  f(l)  so  that  results  a  and  b  follow  immediately. 

2 

Note  that  convergence  for  rotation-reflection  is  o(l/n  )  as  compared 
2 

to  0(1 /n  )  for  rotation  sampling  alone. 


Corollary  10.1  .  If  g  €  C2[0,1]  ,  then  var  Tn  =  0(l/n^)  . 


Proof  .  Expression  (16)  gives 

B2(V)[f'(l)  -  f'(0)] 


*★ 

Tn  =  $  +  4 


T 


+  o(l/n2) 
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from  which  var  Tn  =  0(l/n4)  .  Note  that  the  additional  smoothness 
in  g  considerably  increases  the  convergence  rate. 

Theorem  11  .  If  g  is  piecewise  linear  with  finite  discontinuities, 
then  var  Tn  =  0(l/n^)  . 

Proof  .  See  the  proof  of  Theorem  6  in  the  Appendix  with  f  =  R  . 
Example  ll.l  .  Consider  the  Bernoulli  case  of  Example  6.1  .  Here 
Q(6)  =  (p  -  | (1  -  0)  -  p|)+  _  *  from  which  it  follows  that  2R(d)  = 

(p  -  e)  +  (p  +  0  -  1)  -  p  =  P(e)  .  Therefore,  var  T^  = 

(np  mod  1)(1  -  (np  mod  l))/n  =  var  Tn  .  Also,  R  is  convex  so  that 
0n  retains  its  optimal  property. 

Unbounded  g  . 

Theorem  12  .  For  g(u)  =  ua(l  -  u)br(u)  where  a,b  <  0  ,  g  is  inte- 
grable  and  r  £  C^[0,1]  ,  var  Tn*  =  0(l/n2(1+al)  for  a  <  b  and 
var  Tn  =  0(l/n2(^+b))  for  a  >  b  . 

Proof  .  By  appropriate  use  of  Lemma  5.1  one  can  show  that 


where  £  =  mU  mod  1  .  Using  (17)  leads  to  the  result.  Note  the 
absence  of  any  advantage  in  terms  of  the  ultimate  convergence  rate. 

Theorem  13  .  For  g(u)  =  ua(l-u)bq(u)  In  u  where  a.b  <  0  ,  g  is 
integrable  and  q  CC^[0,1]  , 

a.  var  Tn  =  0((ln  n/n^+a)^)  a  <  b  <  0 

b.  var  T**  =  0(l/n2(1+b))  b<a 

c.  var  Tn  =  0(l/n4) 


a  =  b  =  0  . 
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Proof  .  Results  a  and  b  follow  directly  from  (18)  and  (19)  as 
in  Theorem  8  .  Result  b  arises  as  a  consequence  of  z;(0,x)  =  -^(0,1-x)  . 
The  implication  of  result  c  for  rotation-reflection  sampling  is  that 
the  large  sample  convergence  rate  0(l/n  )  is  achieved  faster  than  in 
the  case  of  rotation  sampling  alone.  A  reexamination  of  the  exponential 
case  in  Example  8.1  illustrates  this  case. 


5.  What  About  Discrete  Event  Simulation? 


In  discrete  event  simulation,  the  sampling  problem  incurred  usually 
is  a  multivariate  one  for  which  it  is  known  that  the  variance  reduction 
properties  for  the  univariate  case  do  not  necessarily  hold.  Although 
this  topic  remains  for  future  research,  at  least  one  important  situation 
that  arises  in  congestion  property  establishes  the  importance  of  studying 
the  univariate  case.  Let  us  return  to  the  single  server  queue  simulation 
of  Section  1  .  Let  B..  =  S . .  -  A.,  for  i  =  l,...,m  waiting  times  on 

I  J  I  J  I  J 

replication  j  =  l,...,n  .  Then  as  the  traffic  intensity  approaches  unity, 

i 


Wi3 


Oj 


k=l 


so  that 


W 

* 


j 


m 


(m  -  i  +  1  )B.. j 


Here  .  becomes  a  sum  of  independent  random  variables  and  if  one  uses 
J 

rotation-reflection  sampling  to  generate  B..  n  for  each  i  one 

A 

can  expect  var  pm  n  to  show  a  convergence  rate  associated  with  the  known 
distribution  of  the  B. .  .  Preliminary  sampling  experiments  confirm  this 

■  J 

result  for  large  traffic  intensities. 
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6.  Conclusions 

The  results  presented  here  extend  those  in  Hammersley  and  Morton 

(1956)  by  showing  the  covariance  structures  induced  by  the  rotation 

and  reflection  (antithetic)  sampling  plans,  deriving  conditions  under 

which  these  sampling  plans  are  optimal  and  by  examining  the  unbounded 

case.  For  the  piecewise  linear  case,  the  results  suggest  that  a  sample 

size  n  can  be  considerably  more  desirable  than  another  n'  although 

n'  >  n  .  The  results  also  show  that  the  benefits  of  reflection  sampling 

arise  principally  for  symmetric  (property  16)  functions.  The  benefit 

★  * 

for  nonsymmetric  unbounded  functions  is  to  speed  the  rate  of  var  T^ 
for  moderate  n  to  the  ultimate  rate  achievable  with  rotation  sampling 
alone.  This  is  clearly  advantageous  when  working  within  a  limited  budget. 
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APPENDIX 

Proof  of  Theorem  2  .  Lemma  A.1  shows  that  (8)  is  a  convex  program¬ 
ming  problem.  Then  we  show  z*  =  ,  . . .  ,  1}  is  a  local  minimum  point. 

Since  a  local  minimum  point  of  a  convex  programming  problem  is  also  a 
global  minimum  point,  z*  is  a  global  minimum  point. 

Lemma  A.1  .  Formulation  (8)  is  a  convex  programming  problem. 

Proof  .  For  every  (i,j)  where  1  <  i  <  n-1  and  i+1  <  j  <  n  , 
define 


Li.j-i)*'1  ‘ 


+  z 


j-1 


where 

n-1 

z  C  Z  =  {(z1,...,zn_1)H  zk  <  1,  zk  >  0  k  =  1 , . . .  ,n-l }  . 

k  - 1 

Here  7L  denotes  the  feasible  region  of  (8)  ,  ,  is  a  linear 

\  *  t  J  *"  *  / 

function  on  ZZ  and 


C(z,  ♦  ...  *  Zj.,)  ■  C(l(1ij.1)({D 

is  a  convex  function  on  Z  .  Since  the  objective  function  w  in  (8) 
is  the  sum  of  convex  functions,  w  is  convex  on  Z  .  Since  the  con¬ 
straints  in  (8)  are  linear,  (8)  is  by  definition  a  convex  programming 
problem. 

Let  m  and  m#  be  positive  integers.  Then,  convexity  gives 

mC(a)  +  m'C(b)  a  (m  +  m')C(m*-  +  ffrb)  (A.1) 


for  any  0  s  a,  b  s  1  . 
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Proof  of  the  Theorem  .  Since  z*  is  an  interior  point  of  TL  ,  there 
exists  an  open  neighborhood  N(z*)  such  that 


z* €  N(z*)  C  TL  . 


Then  z*  is  a  local  minimum  point  if 


w(z*)  <  w(z*  +  y) 


(A. 2) 


for  all  z*  +  y  in  N(z*)  .  The  value  of  the  objective  function  at 
the  perturbed  point  z*  +  y  is  from  (8) 


w(z  +  y )  ~  wC(z-)  )  +  (y i » . . .  »y^_ i ) J 

=  w[l+  yr  ....  J+  yn.-|] 

n-1  n-j  .  i+j-1 

-l  l  C(i  *1  yt)  . 
j=l  i=l  "  t=i  1 


Rearranging  the  terms  in  (A. 3)  ,  we  have 


(A. 3) 


»U*  *  y>  -  tr;  C(i  *  y()  +  C(!tl  +  y,)] 

2 


(A. 4) 


n-2  n  2  n  ?  i+n_1 

+  [IW  +yi+1)  *  C(V+  It„  V^ 


nC(l)  +  nC(|)  +  ...  +  nC(^)  n  odd  (A. 5) 

nC(lj  +  nC(£)  +  ...  +  £  C(l)  n  even  . 
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The  last  equality  of  (A. 5)  follows  from  the  symmetry  of  C  .  That  is. 


C(£)  =  C(^) 


k  =  1 .  ,n-l  . 


Now,  we  show  w(z*  +  y)  >w(z*)  by  proving  that 

n-k  .  i+k-1  k  „  .  n+i-k-1 

l  c(i  +  l  Jt)  *  l  *  l  yt) 

i=l  n  t=i  1  i=l  n  t=i  1 

n-k  .  k  .  . 

a  I  C(£)  +  l  C(^)  =  nC(J)  k  =  1 . n-1  . 

i*l  n  i*l  n  n 


(A. 6) 


i+k-1 


Let  r.  .  =  l  y  Repeatedly  using  (A. 2)  ,  we  have 
1,K  t=i  1 

^n-k 

n-k  v  .  ,  ri,k 


(A. 7) 


The  symmetry  of  C  gives 


zi-i c<1^  *  r^k> "  zi-i e<-  ‘  “ 


l  ri  k 


(A. 8) 


The  last  inequality  of  (A. 8)  also  follows  from  repeated  use  of  (A.l)  . 
Combining  the  results  of  (A. 7)  and  (A. 8)  ,  we  have 

Cc("tr+k)+Cc(^tr'--k) 

n-k  k 

2  <"  -  k>c<!f  *  ) 2  "c<!;>  ■ 


which  proves  the  inequality  (A. 6)  .  From  (A. 6),  (A. 5),  and  (A. 4)  we  have 
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w(z*  +  y)  >  w(z*) 

for  all  z*  +  y  €  N(z*)  ,  which  completes  the  proof. 


Proof  of  Theorem  6  .  The  proof  follows  as  a  consequence  of  Lemma  A. 2  . 

Lemma  A. 2  .  Let  f  be  a  continuous  piecewise  linear  function  on  [0,1] 
with  parameters  s-j ,  S2,  and  c  such  that 

f ( x )  =  [s.,x  +  d-, ) 1 C 0 , c ) ( X )  +  ^S2X  +  d2 ^ 1  [ c ,  1  1^  0  '  c  ^  1  , 

where  I  denotes  the  indicator  function.  Then  the  quantity 


e 

n 


f(0)  +  f(D  +  £ 


1 

/  f(x)  dx 

0 


2 

decreases  as  0(l/n  )  . 

Proof  .  Given  n,  let  [ i n/n ,  (i  +l)/n]  be  the  subinterval  that 
includes  c  .  Then 


(continued  next  page) 
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S2  -  S1 


i  i  . 

Since  —  <  c  <  +  I  ,  we  have 


and 


in  _  nc  -  LncJ  1 


V1 


-  c 


V1  - nc 


1  -  {nc  -  LncJ  ) 


<  1 

n 


Therefore, 


e„<f) 


s2  '  S1 


2n 


1 —  (nc  -  LncJ  ) (1  -  (nc  -  Lncji) 


which  is  zero  or  decreases  as  0(l/n^)  . 

Extension  to  a  general  continuous  piecewise  linear  function 


=  2  1  fc  c  ~  +  d,] 

i=l  Lci’ci+1;  1  1  1 


is  direct. 

To  prove  Theorem  6  one  needs  only  set  f(x)  =  P(x)  for  0  <  x  <  1 
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